require(estimatr)

#### data preparation ####
FtF.data <- read.csv("FtF_data.csv")
nrow(FtF.data)  # number of respondents

online.data <- read.csv("online_data.csv")
nrow(online.data)  # number of respondents

## young-respondent subset of the face-to-face data
FtF.young.data <- subset(FtF.data, age < 45)

#### test of Hypothesis 1 (FtF) ####
## DQ-based estimate
DQ.est.FtF <- mean(FtF.data$DQ, na.rm = TRUE)
DQ.se.FtF <- sqrt(DQ.est.FtF * (1 - DQ.est.FtF) / sum(! is.na(FtF.data$DQ)))
round(DQ.est.FtF, 3)
round(DQ.est.FtF + qnorm(0.025) * DQ.se.FtF, 3)
round(DQ.est.FtF + qnorm(0.975) * DQ.se.FtF, 3)

## list-based estimate
list.result.FtF <- lm_robust(list ~ treatment, data = FtF.data, 
                             fixed_effects = ~ block)
list.est.FtF <- as.numeric(list.result.FtF$coefficients)
list.se.FtF <- as.numeric(list.result.FtF$std.error)
round(list.est.FtF, 3)
round(list.est.FtF + qnorm(0.025) * list.se.FtF, 3)
round(list.est.FtF + qnorm(0.975) * list.se.FtF, 3)

## SDB estimate
SDB.est.FtF <- DQ.est.FtF - list.est.FtF
SDB.se.FtF <- sqrt(list.se.FtF ^ 2 + DQ.se.FtF ^ 2)
round(SDB.est.FtF, 3)
round(SDB.est.FtF + qnorm(0.025) * SDB.se.FtF, 3)
round(SDB.est.FtF + qnorm(0.975) * SDB.se.FtF, 3)
round((1 - pnorm(abs(SDB.est.FtF / SDB.se.FtF))) * 2, 3)

## Pulse Asia-style DQ question (note 16)
round(mean(FtF.data$PulseAsia > 3), 3)
round(mean(FtF.data$PulseAsia == 3), 3)  # middle category

#### test of Hypothesis 1 (online) ####
## DQ-based estimate
DQ.est.online <- mean(online.data$DQ, na.rm = TRUE)
DQ.se.online <- sqrt(DQ.est.online * (1 - DQ.est.online) / 
                             sum(! is.na(online.data$DQ)))
round(DQ.est.online, 3)
round(DQ.est.online + qnorm(0.025) * DQ.se.online, 3)
round(DQ.est.online + qnorm(0.975) * DQ.se.online, 3)

## list-based estimate
list.result.online <- lm_robust(list ~ treatment, data = online.data, 
                                fixed_effects = ~ block)
list.est.online <- as.numeric(list.result.online$coefficients)
list.se.online <- as.numeric(list.result.online$std.error)
round(list.est.online, 3)
round(list.est.online + qnorm(0.025) * list.se.online, 3)
round(list.est.online + qnorm(0.975) * list.se.online, 3)

## SDB estimate
SDB.est.online <- DQ.est.online - list.est.online
SDB.se.online <- sqrt(list.se.online ^ 2 + DQ.se.online ^ 2)
round(SDB.est.online, 3)
round(SDB.est.online + qnorm(0.025) * SDB.se.online, 3)
round(SDB.est.online + qnorm(0.975) * SDB.se.online, 3)
round((1 - pnorm(abs(SDB.est.online / SDB.se.online))) * 2, 3)

## Pulse Asia-style DQ question (note 16)
round(mean(online.data$PulseAsia > 3, na.rm = TRUE), 3)
round(mean(online.data$PulseAsia == 3, na.rm = TRUE), 3)  # middle category

#### test of Hypothesis 2 ####
## SDB estimate using the young-respondent FtF data
DQ.est.FtF.young <- mean(FtF.young.data$DQ, na.rm = TRUE)
DQ.se.FtF.young <- sqrt(DQ.est.FtF.young * (1 - DQ.est.FtF.young) / 
                                sum(! is.na(FtF.young.data$DQ)))

list.result.FtF.young <- lm_robust(list ~ treatment, data = FtF.young.data, 
                                   fixed_effects = ~ block)
list.est.FtF.young <- as.numeric(list.result.FtF.young$coefficients)
list.se.FtF.young <- as.numeric(list.result.FtF.young$std.error)

SDB.est.FtF.young <- DQ.est.FtF.young - list.est.FtF.young
SDB.se.FtF.young <- sqrt(list.se.FtF.young ^ 2 + DQ.se.FtF.young ^ 2)

## difference in SDB between survey modes
SDB.diff.est <- SDB.est.FtF.young - SDB.est.online
SDB.diff.se <- sqrt(list.se.FtF.young ^ 2 + DQ.se.FtF.young ^ 2 + 
                      list.se.online ^ 2 + DQ.se.online ^ 2)
round(SDB.diff.est, 3)
round(SDB.diff.est + qnorm(0.025) * SDB.diff.se, 3)
round(SDB.diff.est + qnorm(0.975) * SDB.diff.se, 3)
round((1 - pnorm(abs(SDB.diff.est / SDB.diff.se))) * 2, 3)

#### Figure 2 ####
cairo_pdf("Figure_2.pdf", width = 4.6, height = 2.6, pointsize = 8)
par(mar = c(4, 0, 3, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-0.7, 1.9), ylim = c(0.5, 4.5), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(seq(0.2, 1, 0.2), seq(1.3, 1.9, 0.2)), col = "gray")
segments(DQ.est.FtF + qnorm(0.025) * DQ.se.FtF, 4.1, 
         DQ.est.FtF + qnorm(0.975) * DQ.se.FtF, 4.1)
points(DQ.est.FtF, 4.1, pch = 15)
text(DQ.est.FtF, 4.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF + qnorm(0.025) * list.se.FtF, 3.9, 
         list.est.FtF + qnorm(0.975) * list.se.FtF, 3.9)
points(list.est.FtF, 3.9, pch = 19)
text(list.est.FtF, 3.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.FtF.young + qnorm(0.025) * DQ.se.FtF.young, 3.1, 
         DQ.est.FtF.young + qnorm(0.975) * DQ.se.FtF.young, 3.1)
points(DQ.est.FtF.young, 3.1, pch = 15)
text(DQ.est.FtF.young, 3.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.FtF.young + qnorm(0.025) * list.se.FtF.young, 2.9, 
         list.est.FtF.young + qnorm(0.975) * list.se.FtF.young, 2.9)
points(list.est.FtF.young, 2.9, pch = 19)
text(list.est.FtF.young, 2.9, "List", pos = 1, cex = 0.8)
segments(DQ.est.online + qnorm(0.025) * DQ.se.online, 2.1, 
         DQ.est.online + qnorm(0.975) * DQ.se.online, 2.1)
points(DQ.est.online, 2.1, pch = 15)
text(DQ.est.online, 2.1, "DQ", pos = 3, cex = 0.8)
segments(list.est.online + qnorm(0.025) * list.se.online, 1.9, 
         list.est.online + qnorm(0.975) * list.se.online, 1.9)
points(list.est.online, 1.9, pch = 19)
text(list.est.online, 1.9, "List", pos = 1, cex = 0.8)
segments(SDB.est.FtF + qnorm(0.025) * SDB.se.FtF + 1.3, 4, 
         SDB.est.FtF + qnorm(0.975) * SDB.se.FtF + 1.3, 4)
points(SDB.est.FtF + 1.3, 4, pch = 4)
text(SDB.est.FtF + 1.3, 4, sprintf("%5.3f", SDB.est.FtF), 
     cex = 0.8, pos = 3)
segments(SDB.est.FtF.young + qnorm(0.025) * SDB.se.FtF.young + 1.3, 3, 
         SDB.est.FtF.young + qnorm(0.975) * SDB.se.FtF.young + 1.3, 3)
points(SDB.est.FtF.young + 1.3, 3, pch = 4)
text(SDB.est.FtF.young + 1.3, 3, sprintf("%5.3f", SDB.est.FtF.young), 
     cex = 0.8, pos = 3)
segments(SDB.est.online + qnorm(0.025) * SDB.se.online + 1.3, 2, 
         SDB.est.online + qnorm(0.975) * SDB.se.online + 1.3, 2)
points(SDB.est.online + 1.3, 2, pch = 4)
text(SDB.est.online + 1.3, 2, sprintf("%5.3f", SDB.est.online), 
     cex = 0.8, pos = 3)
segments(SDB.diff.est + qnorm(0.025) * SDB.diff.se + 1.3, 1, 
         SDB.diff.est + qnorm(0.975) * SDB.diff.se + 1.3, 1)
points(SDB.diff.est + 1.3, 1, pch = 4)
text(SDB.diff.est + 1.3, 1, sprintf("%5.3f", SDB.diff.est), 
     cex = 0.8, pos = 3)
text(-0.7, 4, "Face-to-face survey\n(entire)", pos = 4)
text(-0.7, 3, "Face-to-face survey\n(18\u201344)", pos = 4)
text(-0.7, 2, "Online survey", pos = 4)
text(-0.7, 1, "Difference in differences\nb/w survey modes", pos = 4)
axis(1, seq(0.2, 1, 0.2), lwd = 0.5)
axis(1, seq(1.3, 1.9, 0.2), labels = c("0.0", "0.2", "0.4", "0.6"), lwd = 0.5)
mtext("Percent", side = 1, line = 2.5, at = 0.6)
mtext("Percentage points", side = 1, line = 2.5, at = 1.6)
mtext("Approval rate\n", line = 0.5, at = 0.6, font = 2, cex = 1.1)
mtext("Social desirability\nbias", line = 0.5, at = 1.6, font = 2, cex = 1.1)
dev.off()